Libraries

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.0
## ✓ tidyr   1.1.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## Warning: package 'ggplot2' was built under R version 3.6.2
## Warning: package 'tibble' was built under R version 3.6.2
## Warning: package 'tidyr' was built under R version 3.6.2
## Warning: package 'purrr' was built under R version 3.6.2
## Warning: package 'dplyr' was built under R version 3.6.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following object is masked from 'package:purrr':
## 
##     transpose
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
## Loading required package: nlme
## Warning: package 'nlme' was built under R version 3.6.2
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
## Loading required package: emmeans
## Warning: package 'emmeans' was built under R version 3.6.2
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## Warning: package 'multcomp' was built under R version 3.6.2
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 3.6.2
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.6.2
## Loading required package: TH.data
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 3.6.2
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Functions

Loading in the data

Demographics

Unresitricted data

Graph data

Individual data

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   group = col_character(),
##   measure = col_character()
## )
## See spec(...) for full column specifications.
## [1] 1029  103
## Parsed with column specification:
## cols(
##   sub = col_double(),
##   FC = col_double(),
##   group = col_character()
## )
## Parsed with column specification:
## cols(
##   sub = col_double(),
##   FC = col_double(),
##   group = col_character()
## )
## Parsed with column specification:
## cols(
##   sub = col_double(),
##   FC = col_double(),
##   group = col_character()
## )
## [1] 343   3
## Joining by: sub, group

joining individual sub data

## Joining by: sub, group
## [1] 1029  118
## Joining by: sub
## [1] 1029  120

Cleanup, center, ect

## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## 'data.frame':    1029 obs. of  15 variables:
##  $ IC_98            : num  0.00412 0.00362 0.00481 0.0027 0.00449 ...
##  $ IC_99            : num  0.00377 0.00454 0.00343 0.00579 0.00428 ...
##  $ group            : Factor w/ 3 levels "no","ob","ov": 1 1 1 1 1 1 1 1 1 1 ...
##  $ measure          : Factor w/ 3 levels "centrality","cluster",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ FC               : num  0.134 0.123 0.151 0.16 0.105 ...
##  $ Age_in_Yrs       : int  22 29 35 27 26 34 33 25 33 31 ...
##  $ ZygosityGT       : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
##  $ Race             : Factor w/ 6 levels "Am. Indian/Alaskan Nat.",..: 6 3 6 6 6 6 6 3 2 6 ...
##  $ Ethnicity        : Factor w/ 3 levels "Hispanic/Latino",..: 2 2 1 2 2 2 2 2 2 2 ...
##  $ BMI              : num  21.1 22.2 20.4 22.6 24.5 ...
##  $ HbA1C            : num  NA 5.4 5.3 5.4 4.7 4.7 5.5 NA NA NA ...
##  $ Hypothyroidism   : int  NA NA 0 NA NA 0 0 0 0 NA ...
##  $ Hyperthyroidism  : int  NA NA 0 NA NA 0 0 0 0 NA ...
##  $ OtherEndocrn_Prob: int  NA NA 0 NA NA 0 0 0 0 NA ...
##  $ Mother_ID        : Factor w/ 339 levels "50263","50371",..: 140 186 141 231 60 278 236 192 262 15 ...
##  [1] ""                    "Amygdala"            "Auditory"           
##  [4] "BrainStem"           "Caudate"             "Cerebellum"         
##  [7] "CinguloOperc"        "Default"             "DorsalAttn"         
## [10] "FrontoParietal"      "MedialParietal"      "None"               
## [13] "Putamen"             "Smhand"              "Smmouth"            
## [16] "Thalamus"            "VentralAttn"         "VentralDiencephalon"
## [19] "Visual"
## The following `from` values were not present in `x`: cerebellum

Contrasts for demographics

## [1] "no" "ov" "ob"
## [1] "White"                                 
## [2] "Black or African American"             
## [3] "Asian"                                 
## [4] "Multi-racial"                          
## [5] "Unknown or not reported and Hispanic"  
## [6] "Multi-racial and Hispanic"             
## [7] "Black or African American and Hispanic"
## [8] "White and Hispanic"

Covariates

Age

##  contrast estimate    SE  df t.ratio p.value
##  noVov       0.134 0.488 340  0.275  0.9898 
##  noVob      -1.019 0.542 340 -1.880  0.1719 
##  ovVob      -1.153 0.580 340 -1.990  0.1356 
## 
## P value adjustment: sidak method for 3 tests

## 
##  Descriptive statistics by group 
## group: no
##    vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 156 28.33 4.11   27.5   28.23 5.19  22  36    14  0.2    -1.27 0.33
## ------------------------------------------------------------ 
## group: ov
##    vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 109 28.19 3.72     28   28.16 4.45  22  35    13 0.02    -1.06 0.36
## ------------------------------------------------------------ 
## group: ob
##    vars  n  mean   sd median trimmed  mad min max range  skew kurtosis   se
## X1    1 78 29.35 3.75     30   29.48 4.45  22  35    13 -0.23     -1.1 0.43

BMI

##  contrast estimate    SE  df t.ratio p.value
##  noVov       -5.01 0.280 340 -17.874 <.0001 
##  noVob      -12.03 0.311 340 -38.653 <.0001 
##  ovVob       -7.02 0.333 340 -21.097 <.0001 
## 
## P value adjustment: sidak method for 3 tests

## 
##  Descriptive statistics by group 
## group: no
##    vars   n mean   sd median trimmed  mad   min   max range  skew kurtosis   se
## X1    1 156 22.2 1.69  22.38   22.25 1.93 18.11 24.95  6.84 -0.28    -0.86 0.14
## ------------------------------------------------------------ 
## group: ov
##    vars   n mean   sd median trimmed  mad min   max range skew kurtosis   se
## X1    1 109 27.2 1.56  27.06   27.15 2.05  25 29.98  4.98 0.21     -1.3 0.15
## ------------------------------------------------------------ 
## group: ob
##    vars  n  mean   sd median trimmed  mad  min   max range skew kurtosis   se
## X1    1 78 34.23 3.62  33.39   33.76 2.79 30.1 43.94 13.84 1.08     0.29 0.41

#HBA1c

##  contrast estimate     SE  df t.ratio p.value
##  noVov     -0.0178 0.0645 224 -0.276  0.9897 
##  noVob     -0.1779 0.0745 224 -2.388  0.0524 
##  ovVob     -0.1601 0.0786 224 -2.037  0.1230 
## 
## P value adjustment: sidak method for 3 tests
## Warning: Removed 116 rows containing non-finite values (stat_ydensity).
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to unique
## 'x' values

## 
##  Descriptive statistics by group 
## group: no
##    vars   n mean   sd median trimmed mad  min max range  skew kurtosis   se
## X1    1 103 5.21 0.36    5.2    5.23 0.3 3.12   6  2.88 -1.83     9.27 0.04
## ------------------------------------------------------------ 
## group: ov
##    vars  n mean   sd median trimmed mad min max range skew kurtosis   se
## X1    1 76 5.23 0.24    5.2    5.23 0.3 4.8 5.9   1.1 0.17    -0.29 0.03
## ------------------------------------------------------------ 
## group: ob
##    vars  n mean  sd median trimmed mad min max range skew kurtosis  se
## X1    1 48 5.39 0.7    5.3    5.33 0.3 4.5 9.6   5.1 4.34    23.93 0.1

sex

## 
##  Pearson's Chi-squared test
## 
## data:  table(dat$group, dat$Gender)
## X-squared = 9.4234, df = 2, p-value = 0.008989
## $F
## no ov ob 
## 90 42 39 
## 
## $M
## no ov ob 
## 66 67 39

Race

## Warning in chisq.test(table(dat$group, dat$race_eth)): Chi-squared approximation
## may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(dat$group, dat$race_eth)
## X-squared = 22.324, df = 14, p-value = 0.07219
## $White
## no ov ob 
## 94 72 51 
## 
## $`Black or African American`
## no ov ob 
## 22 17 14 
## 
## $Asian
## no ov ob 
## 19  6  0 
## 
## $`Multi-racial`
## no ov ob 
##  4  3  1 
## 
## $`Unknown or not reported and Hispanic`
## no ov ob 
##  0  2  4 
## 
## $`Multi-racial and Hispanic`
## no ov ob 
##  2  1  1 
## 
## $`Black or African American and Hispanic`
## no ov ob 
##  1  0  0 
## 
## $`White and Hispanic`
## no ov ob 
## 14  7  7

Models

fit_function<-function(y,data){
  fit<-lm(y~group+Age_c+HBA1c_c+Gender+race_eth+FC, data=data)
  test.emm <- emmeans(fit, "group")
  contest<-emmeans(test.emm, pairwise ~ group)
  test.emm.s <- update(contest$contrasts, infer = c(TRUE, TRUE))
  return(summary(test.emm.s))
}

dfer<-function(X,Y){
  df<-data.frame(X,Y, stringsAsFactors=FALSE)
  return(df)
}

clustering

cluster<-subset(data, data$measure == "cluster")

clus_vals<-apply(cluster[2:101], FUN=fit_function, MARGIN = 2 ,data=cluster)
indx<-names(clus_vals)

l <- list()
for (i in seq(1, length(clus_vals), by= 1)){
  print(paste("start", i))
  df <- dfer(indx[i],clus_vals[i])
  colnames(df) <- c("IC", "contrast", "estimate","SE","df", "lower.CL", "upper.CL", "t.ratio", "p.value")
  l[[i]] <-df
}
## [1] "start 1"
## [1] "start 2"
## [1] "start 3"
## [1] "start 4"
## [1] "start 5"
## [1] "start 6"
## [1] "start 7"
## [1] "start 8"
## [1] "start 9"
## [1] "start 10"
## [1] "start 11"
## [1] "start 12"
## [1] "start 13"
## [1] "start 14"
## [1] "start 15"
## [1] "start 16"
## [1] "start 17"
## [1] "start 18"
## [1] "start 19"
## [1] "start 20"
## [1] "start 21"
## [1] "start 22"
## [1] "start 23"
## [1] "start 24"
## [1] "start 25"
## [1] "start 26"
## [1] "start 27"
## [1] "start 28"
## [1] "start 29"
## [1] "start 30"
## [1] "start 31"
## [1] "start 32"
## [1] "start 33"
## [1] "start 34"
## [1] "start 35"
## [1] "start 36"
## [1] "start 37"
## [1] "start 38"
## [1] "start 39"
## [1] "start 40"
## [1] "start 41"
## [1] "start 42"
## [1] "start 43"
## [1] "start 44"
## [1] "start 45"
## [1] "start 46"
## [1] "start 47"
## [1] "start 48"
## [1] "start 49"
## [1] "start 50"
## [1] "start 51"
## [1] "start 52"
## [1] "start 53"
## [1] "start 54"
## [1] "start 55"
## [1] "start 56"
## [1] "start 57"
## [1] "start 58"
## [1] "start 59"
## [1] "start 60"
## [1] "start 61"
## [1] "start 62"
## [1] "start 63"
## [1] "start 64"
## [1] "start 65"
## [1] "start 66"
## [1] "start 67"
## [1] "start 68"
## [1] "start 69"
## [1] "start 70"
## [1] "start 71"
## [1] "start 72"
## [1] "start 73"
## [1] "start 74"
## [1] "start 75"
## [1] "start 76"
## [1] "start 77"
## [1] "start 78"
## [1] "start 79"
## [1] "start 80"
## [1] "start 81"
## [1] "start 82"
## [1] "start 83"
## [1] "start 84"
## [1] "start 85"
## [1] "start 86"
## [1] "start 87"
## [1] "start 88"
## [1] "start 89"
## [1] "start 90"
## [1] "start 91"
## [1] "start 92"
## [1] "start 93"
## [1] "start 94"
## [1] "start 95"
## [1] "start 96"
## [1] "start 97"
## [1] "start 98"
## [1] "start 99"
## [1] "start 100"
df_cluster<-bind_rows(l)
df_cluster$pFDR<-p.adjust(df_cluster$p.value, method = "BH")

df_cluster[df_cluster$p.value <= 0.05, ] 
violin(cluster,cluster$group,cluster$IC_16)

violin(cluster,cluster$group,cluster$IC_23)

centrality

centrality<-subset(data, data$measure == "centrality")

cen_vals<-apply(centrality[2:101], FUN=fit_function, MARGIN = 2 ,data=centrality)
indx<-names(cen_vals)

l <- list()
for (i in seq(1, length(cen_vals), by= 1)){
  print(paste("start", i))
  df <- dfer(indx[i],cen_vals[i])
  colnames(df) <- c("IC", "contrast", "estimate","SE","df", "lower.CL", "upper.CL", "t.ratio", "p.value")
  l[[i]] <-df
}
## [1] "start 1"
## [1] "start 2"
## [1] "start 3"
## [1] "start 4"
## [1] "start 5"
## [1] "start 6"
## [1] "start 7"
## [1] "start 8"
## [1] "start 9"
## [1] "start 10"
## [1] "start 11"
## [1] "start 12"
## [1] "start 13"
## [1] "start 14"
## [1] "start 15"
## [1] "start 16"
## [1] "start 17"
## [1] "start 18"
## [1] "start 19"
## [1] "start 20"
## [1] "start 21"
## [1] "start 22"
## [1] "start 23"
## [1] "start 24"
## [1] "start 25"
## [1] "start 26"
## [1] "start 27"
## [1] "start 28"
## [1] "start 29"
## [1] "start 30"
## [1] "start 31"
## [1] "start 32"
## [1] "start 33"
## [1] "start 34"
## [1] "start 35"
## [1] "start 36"
## [1] "start 37"
## [1] "start 38"
## [1] "start 39"
## [1] "start 40"
## [1] "start 41"
## [1] "start 42"
## [1] "start 43"
## [1] "start 44"
## [1] "start 45"
## [1] "start 46"
## [1] "start 47"
## [1] "start 48"
## [1] "start 49"
## [1] "start 50"
## [1] "start 51"
## [1] "start 52"
## [1] "start 53"
## [1] "start 54"
## [1] "start 55"
## [1] "start 56"
## [1] "start 57"
## [1] "start 58"
## [1] "start 59"
## [1] "start 60"
## [1] "start 61"
## [1] "start 62"
## [1] "start 63"
## [1] "start 64"
## [1] "start 65"
## [1] "start 66"
## [1] "start 67"
## [1] "start 68"
## [1] "start 69"
## [1] "start 70"
## [1] "start 71"
## [1] "start 72"
## [1] "start 73"
## [1] "start 74"
## [1] "start 75"
## [1] "start 76"
## [1] "start 77"
## [1] "start 78"
## [1] "start 79"
## [1] "start 80"
## [1] "start 81"
## [1] "start 82"
## [1] "start 83"
## [1] "start 84"
## [1] "start 85"
## [1] "start 86"
## [1] "start 87"
## [1] "start 88"
## [1] "start 89"
## [1] "start 90"
## [1] "start 91"
## [1] "start 92"
## [1] "start 93"
## [1] "start 94"
## [1] "start 95"
## [1] "start 96"
## [1] "start 97"
## [1] "start 98"
## [1] "start 99"
## [1] "start 100"
df_centrality<-bind_rows(l)
df_centrality$pFDR<-p.adjust(df_centrality$p.value, method = "BH")

df_centrality[df_centrality$p.value <= 0.05, ] 
violin(centrality,centrality$group,centrality$IC_16)

violin(centrality,centrality$group,centrality$IC_11)

participation coeffiecent

PC<-subset(data, data$measure == "PC")

PC_vals<-apply(PC[2:101], FUN=fit_function, MARGIN = 2 ,data=PC)
indx<-names(PC_vals)

l <- list()
for (i in seq(1, length(PC_vals), by= 1)){
  # print(paste("start", i))
  df <- dfer(indx[i],PC_vals[i])
  colnames(df) <- c("IC", "contrast", "estimate","SE","df", "lower.CL", "upper.CL", "t.ratio", "p.value")
  l[[i]] <-df
}

df_PC<-bind_rows(l)
df_PC$pFDR<-p.adjust(df_PC$p.value, method = "BH")

df_PC[df_PC$p.value <= 0.05, ] 
chisq.test(table(graph_df$group, graph_df$node_type))
## Warning in chisq.test(table(graph_df$group, graph_df$node_type)): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(graph_df$group, graph_df$node_type)
## X-squared = 4.7689, df = 4, p-value = 0.3118
tapply(graph_df$group, graph_df$node_type, summary)
## $connector
##     normal overweight      obese 
##         34         32         37 
## 
## $kinless
##     normal overweight      obese 
##          0          0          2 
## 
## $peripheral
##     normal overweight      obese 
##         66         68         61
balloonplot(table(graph_df$group, graph_df$node_type), main ="Node Frequencies", xlab ="", ylab="",
            label = TRUE, show.margins = FALSE, colsrt = 90,text.size=1.5 ,dotsize = 5, colmar=5, dotcolor=cbPalette)

no_<-subset(graph_df, graph_df$group == "normal")
pNo<-balloonplot(table(no_$modules, no_$area), main ="Module Frequencies in areas normal", xlab ="", ylab="",
            label = TRUE, show.margins = FALSE, colsrt = 90, text.size=1 ,dotsize = 5, colmar=5, dotcolor=cbPalette[1])

balloonplot(table(no_$module, no_$IC), main ="Module Frequencies by IC", xlab ="", ylab="",
            label = TRUE, show.margins = FALSE, colsrt = 90, text.size=1 ,dotsize = 5, colmar=5, dotcolor=cbPalette[1])

ov_<-subset(graph_df, graph_df$group == "overweight")
pOv<-balloonplot(table(ov_$modules, ov_$area), main ="Module Frequencies in areas overweight", xlab ="", ylab="",
            label = TRUE, show.margins = FALSE, colsrt = 90, text.size=1 ,dotsize = 5, colmar=5, dotcolor=cbPalette[2])

balloonplot(table(ov_$module, ov_$IC), main ="Module Frequencies by IC", xlab ="", ylab="",
            label = TRUE, show.margins = FALSE, colsrt = 90, text.size=1 ,dotsize = 5, colmar=5, dotcolor=cbPalette[2])

ob_<-subset(graph_df, graph_df$group == "obese")
pOb<-balloonplot(table(ob_$modules, ob_$area), main ="Module Frequencies in areas obese", xlab ="", ylab="",
            label = TRUE, show.margins = FALSE, colsrt = 90, text.size=1 ,dotsize = 5, colmar=5, dotcolor=cbPalette[3])

balloonplot(table(ob_$module, ob_$IC), main ="Module Frequencies by IC", xlab ="", ylab="",
            label = TRUE, show.margins = FALSE, colsrt = 90, text.size=1 ,dotsize = 5, colmar=5, dotcolor=cbPalette[3])

Module summaries Average weight 0 - Visual 1 - Fronto-medial parietal and cerebellum and DMN 2 - Amygdala, auditory, hand mouth, and visual 3 - Cerebellum, CO, dorsal attention 4 - Cerebellum 5 - cerebellum, caudate, putamen, thalamus, VDia 6 - brain stem 7 - brain stem, thalamus 8 - brain stem

High weight 0 - Visual 1 - Amygdala, cerebellum, DMN, fronto-parietal 2 - Visual, hand mouth, dorsal attn 3 - Cerebellum, CO, Fronto-parietal 4 - Cerebellum 5 - caudate, putamen, thalamus, VDia 6 - Brain stem 7 - brain stem

Very high weight 0 - Visual 1 - Cerebellum, DMN, fronto-parietal 2 - Auditory, brainstem, CO, DAttn, Visual, hand mouth 3 - Cerebelllum, CO 4 - Amygdala, cerebellum, DMN 5 - Fronto-parietal 6 - Caudate, putamen, thalamus 7 - brainstem 8 - brain stem, VDia

graph_df$modules<-as.factor(graph_df$modules)
# head(graph_df)

mod_function<-function(y,data){
  fit<-lm(y~group, data=data)
  test.emm <- emmeans(fit, "group")
  contest<-emmeans(test.emm, pairwise ~ group)
  test.emm.s <- update(contest$contrasts, infer = c(TRUE, TRUE))
  return(summary(test.emm.s))
}

mod_function(graph_df$centrality,graph_df)
mod_function(graph_df$clustering,graph_df)
mod_function(graph_df$PC,graph_df)
mod_function(abs(graph_df$zDegree),graph_df)
violin(graph_df, graph_df$group, abs(graph_df$zDegree))

violin(graph_df, graph_df$group, graph_df$centrality)

violin(graph_df, graph_df$group, graph_df$clustering)

violin(graph_df, graph_df$group, graph_df$PC)

# levels(graph_df$area)
graph_df$modules<-as.factor(graph_df$modules)
warp <- transform(graph_df, mods = interaction(group, area))
zD.gls <- gls(zDegree ~ mods, data = warp)
zD.emm <- emmeans(zD.gls, "mods")
zD.fac <- update(zD.emm, levels = list(
                group = c("normal", "overweight", "obese"), 
                area = c("", "Amygdala", "Auditory", "BrainStem", "Caudate", "Cerebellum", "CinguloOperc",
                         "Default", "DorsalAttn", "FrontoParietal", "MedialParietal", "None", "Putamen",
                         "Smhand", "Smmouth", "Thalamus", "VentralAttn", "VentralDiencephalon", "Visual")))

contrast(zD.fac, "pairwise", by = "area", adjust = "fdr")
## area = :
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.45599 0.880 243 -0.518  0.7774 
##  normal - obese       0.24919 0.880 243  0.283  0.7774 
##  overweight - obese   0.70519 0.880 243  0.801  0.7774 
## 
## area = Amygdala:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.84578 1.245 243 -0.679  0.7464 
##  normal - obese       0.03058 1.245 243  0.025  0.9804 
##  overweight - obese   0.87636 1.245 243  0.704  0.7464 
## 
## area = Auditory:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.12536 1.245 243 -0.101  0.9199 
##  normal - obese       0.24051 1.245 243  0.193  0.9199 
##  overweight - obese   0.36586 1.245 243  0.294  0.9199 
## 
## area = BrainStem:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.03087 0.440 243 -0.070  0.9442 
##  normal - obese       1.37744 0.440 243  3.129  0.0030 
##  overweight - obese   1.40830 0.440 243  3.199  0.0030 
## 
## area = Caudate:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight  0.48830 0.880 243  0.555  0.5797 
##  normal - obese      -0.77074 0.880 243 -0.875  0.5733 
##  overweight - obese  -1.25904 0.880 243 -1.430  0.4620 
## 
## area = Cerebellum:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.65542 0.220 243 -2.978  0.0032 
##  normal - obese       1.05901 0.220 243  4.811  <.0001 
##  overweight - obese   1.71442 0.220 243  7.789  <.0001 
## 
## area = CinguloOperc:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -1.12981 0.719 243 -1.572  0.3520 
##  normal - obese      -0.58586 0.719 243 -0.815  0.4500 
##  overweight - obese   0.54395 0.719 243  0.757  0.4500 
## 
## area = Default:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.68551 0.471 243 -1.457  0.2197 
##  normal - obese       0.09304 0.471 243  0.198  0.8434 
##  overweight - obese   0.77855 0.471 243  1.654  0.2197 
## 
## area = DorsalAttn:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.93319 0.623 243 -1.499  0.2028 
##  normal - obese      -0.99171 0.623 243 -1.593  0.2028 
##  overweight - obese  -0.05852 0.623 243 -0.094  0.9252 
## 
## area = FrontoParietal:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -1.05449 0.508 243 -2.075  0.0586 
##  normal - obese      -1.12517 0.508 243 -2.214  0.0586 
##  overweight - obese  -0.07068 0.508 243 -0.139  0.8895 
## 
## area = MedialParietal:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.40920 1.245 243 -0.329  0.7427 
##  normal - obese      -1.45150 1.245 243 -1.166  0.6050 
##  overweight - obese  -1.04230 1.245 243 -0.837  0.6050 
## 
## area = None:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.67461 1.245 243 -0.542  0.8827 
##  normal - obese       0.03058 1.245 243  0.025  0.9804 
##  overweight - obese   0.70519 1.245 243  0.566  0.8827 
## 
## area = Putamen:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight  0.48830 0.719 243  0.679  0.4976 
##  normal - obese      -0.77074 0.719 243 -1.072  0.4271 
##  overweight - obese  -1.25904 0.719 243 -1.751  0.2434 
## 
## area = Smhand:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.12536 0.623 243 -0.201  0.8406 
##  normal - obese      -0.26264 0.623 243 -0.422  0.8406 
##  overweight - obese  -0.13729 0.623 243 -0.221  0.8406 
## 
## area = Smmouth:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.12536 0.880 243 -0.142  0.8869 
##  normal - obese       0.24051 0.880 243  0.273  0.8869 
##  overweight - obese   0.36586 0.880 243  0.416  0.8869 
## 
## area = Thalamus:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight  0.18147 0.623 243  0.292  0.7709 
##  normal - obese      -0.26511 0.623 243 -0.426  0.7709 
##  overweight - obese  -0.44659 0.623 243 -0.717  0.7709 
## 
## area = VentralAttn:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -1.28828 0.880 243 -1.463  0.4341 
##  normal - obese      -0.36761 0.880 243 -0.418  0.6766 
##  overweight - obese   0.92067 0.880 243  1.046  0.4451 
## 
## area = VentralDiencephalon:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight -0.02026 0.880 243 -0.023  0.9817 
##  normal - obese       1.70713 0.880 243  1.939  0.0805 
##  overweight - obese   1.72739 0.880 243  1.962  0.0805 
## 
## area = Visual:
##  contrast            estimate    SE  df t.ratio p.value
##  normal - overweight  0.20103 0.321 243  0.625  0.7985 
##  normal - obese       0.20218 0.321 243  0.629  0.7985 
##  overweight - obese   0.00115 0.321 243  0.004  0.9971 
## 
## Degrees-of-freedom method: df.error 
## P value adjustment: fdr method for 3 tests
cer_z<-subset(warp, warp$area == "Cerebellum")
violin(cer_z,cer_z$group,cer_z$zDegree)

VD_z<-subset(warp, warp$area == "VentralDiencephalon")
describeBy(VD_z$zDegree, VD_z$group)
## 
##  Descriptive statistics by group 
## group: normal
##    vars n mean sd median trimmed mad  min  max range skew kurtosis se
## X1    1 2 1.45  0   1.45    1.45   0 1.45 1.45     0  NaN      NaN  0
## ------------------------------------------------------------ 
## group: overweight
##    vars n mean sd median trimmed mad  min  max range skew kurtosis se
## X1    1 2 1.47  0   1.47    1.47   0 1.47 1.47     0  NaN      NaN  0
## ------------------------------------------------------------ 
## group: obese
##    vars n  mean sd median trimmed mad   min   max range skew kurtosis se
## X1    1 2 -0.26  0  -0.26   -0.26   0 -0.26 -0.26     0  NaN      NaN  0

Cerebellum Ventral Diencephalon

PC.gls <- gls(PC ~ mods, data = warp)
PC.emm <- emmeans(PC.gls, "mods")
PC.fac <- update(PC.emm, levels = list(
                group = c("normal", "overweight", "obese"), 
                area = c("","Amygdala", "Auditory", "BrainStem", "Caudate", "Cerebellum", "CinguloOperc",
                         "Default", "DorsalAttn", "FrontoParietal", "MedialParietal", "None", "Putamen",
                         "Smhand", "Smmouth", "Thalamus", "VentralAttn", "VentralDiencephalon", "Visual")))

contrast(PC.fac, "pairwise", by = "area", adjust = "fdr")
## area = :
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight  0.14157 0.1195 243  1.184  0.3561 
##  normal - obese      -0.04981 0.1195 243 -0.417  0.6772 
##  overweight - obese  -0.19138 0.1195 243 -1.601  0.3319 
## 
## area = Amygdala:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.17106 0.1690 243 -1.012  0.4688 
##  normal - obese       0.05474 0.1690 243  0.324  0.7463 
##  overweight - obese   0.22580 0.1690 243  1.336  0.4688 
## 
## area = Auditory:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.06154 0.1690 243 -0.364  0.7161 
##  normal - obese      -0.30120 0.1690 243 -1.782  0.2280 
##  overweight - obese  -0.23966 0.1690 243 -1.418  0.2363 
## 
## area = BrainStem:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight  0.00549 0.0598 243  0.092  0.9269 
##  normal - obese       0.05289 0.0598 243  0.885  0.6427 
##  overweight - obese   0.04740 0.0598 243  0.793  0.6427 
## 
## area = Caudate:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.07927 0.1195 243 -0.663  0.8583 
##  normal - obese      -0.05791 0.1195 243 -0.485  0.8583 
##  overweight - obese   0.02136 0.1195 243  0.179  0.8583 
## 
## area = Cerebellum:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight  0.04365 0.0299 243  1.461  0.2180 
##  normal - obese      -0.00996 0.0299 243 -0.333  0.7391 
##  overweight - obese  -0.05361 0.0299 243 -1.794  0.2180 
## 
## area = CinguloOperc:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.07612 0.0976 243 -0.780  0.8183 
##  normal - obese      -0.05367 0.0976 243 -0.550  0.8183 
##  overweight - obese   0.02244 0.0976 243  0.230  0.8183 
## 
## area = Default:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight  0.07064 0.0639 243  1.106  0.4744 
##  normal - obese       0.06416 0.0639 243  1.004  0.4744 
##  overweight - obese  -0.00648 0.0639 243 -0.102  0.9192 
## 
## area = DorsalAttn:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.13165 0.0845 243 -1.558  0.1809 
##  normal - obese      -0.20977 0.0845 243 -2.482  0.0412 
##  overweight - obese  -0.07813 0.0845 243 -0.924  0.3562 
## 
## area = FrontoParietal:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.09754 0.0690 243 -1.413  0.1588 
##  normal - obese      -0.20679 0.0690 243 -2.997  0.0090 
##  overweight - obese  -0.10925 0.0690 243 -1.583  0.1588 
## 
## area = MedialParietal:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.07480 0.1690 243 -0.442  0.8968 
##  normal - obese      -0.09675 0.1690 243 -0.572  0.8968 
##  overweight - obese  -0.02196 0.1690 243 -0.130  0.8968 
## 
## area = None:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.06554 0.1690 243 -0.388  0.7894 
##  normal - obese       0.04520 0.1690 243  0.267  0.7894 
##  overweight - obese   0.11074 0.1690 243  0.655  0.7894 
## 
## area = Putamen:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.08323 0.0976 243 -0.853  0.7410 
##  normal - obese      -0.03230 0.0976 243 -0.331  0.7410 
##  overweight - obese   0.05094 0.0976 243  0.522  0.7410 
## 
## area = Smhand:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.07013 0.0845 243 -0.830  0.4075 
##  normal - obese      -0.26769 0.0845 243 -3.167  0.0052 
##  overweight - obese  -0.19756 0.0845 243 -2.338  0.0303 
## 
## area = Smmouth:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.01646 0.1195 243 -0.138  0.8906 
##  normal - obese      -0.23472 0.1195 243 -1.964  0.1036 
##  overweight - obese  -0.21826 0.1195 243 -1.826  0.1036 
## 
## area = Thalamus:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight -0.05673 0.0845 243 -0.671  0.5059 
##  normal - obese      -0.11304 0.0845 243 -1.337  0.5059 
##  overweight - obese  -0.05631 0.0845 243 -0.666  0.5059 
## 
## area = VentralAttn:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight  0.04643 0.1195 243  0.388  0.9081 
##  normal - obese       0.01381 0.1195 243  0.116  0.9081 
##  overweight - obese  -0.03262 0.1195 243 -0.273  0.9081 
## 
## area = VentralDiencephalon:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight  0.01901 0.1195 243  0.159  0.8738 
##  normal - obese       0.15633 0.1195 243  1.308  0.3776 
##  overweight - obese   0.13732 0.1195 243  1.149  0.3776 
## 
## area = Visual:
##  contrast            estimate     SE  df t.ratio p.value
##  normal - overweight  0.08225 0.0436 243  1.885  0.0607 
##  normal - obese      -0.13303 0.0436 243 -3.048  0.0038 
##  overweight - obese  -0.21528 0.0436 243 -4.933  <.0001 
## 
## Degrees-of-freedom method: df.error 
## P value adjustment: fdr method for 3 tests
DMN_PC<-subset(warp, warp$area == "Default")
violin(DMN_PC, DMN_PC$group, DMN_PC$PC)

FP_PPC<-subset(warp, warp$area == "FrontoParietal")
violin(FP_PPC, FP_PPC$group, FP_PPC$PC)

vis_PPC<-subset(warp, warp$area == "Visual")
violin(vis_PPC, vis_PPC$group, vis_PPC$PC)

# head(warp)
clus.gls <- gls(clustering ~ mods, data = warp)
clus.emm <- emmeans(clus.gls, "mods")
clus.fac <- update(clus.emm, levels = list(
                group = c("normal", "overweight", "obese"), 
                area = c("","Amygdala", "Auditory", "BrainStem", "Caudate", "Cerebellum", "CinguloOperc",
                         "Default", "DorsalAttn", "FrontoParietal", "MedialParietal", "None", "Putamen",
                         "Smhand", "Smmouth", "Thalamus", "VentralAttn", "VentralDiencephalon", "Visual")))

contrast(clus.fac, "pairwise", by = "area", adjust = "fdr")
## area = :
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.005474 0.0575 243  0.095  0.9920 
##  normal - obese       0.004899 0.0575 243  0.085  0.9920 
##  overweight - obese  -0.000575 0.0575 243 -0.010  0.9920 
## 
## area = Amygdala:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight -0.028490 0.0814 243 -0.350  0.8267 
##  normal - obese      -0.046325 0.0814 243 -0.569  0.8267 
##  overweight - obese  -0.017835 0.0814 243 -0.219  0.8267 
## 
## area = Auditory:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.028521 0.0814 243  0.350  0.8836 
##  normal - obese       0.016597 0.0814 243  0.204  0.8836 
##  overweight - obese  -0.011924 0.0814 243 -0.147  0.8836 
## 
## area = BrainStem:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.013192 0.0288 243  0.459  0.8213 
##  normal - obese       0.019700 0.0288 243  0.685  0.8213 
##  overweight - obese   0.006507 0.0288 243  0.226  0.8213 
## 
## area = Caudate:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.025994 0.0575 243  0.452  0.9778 
##  normal - obese       0.026241 0.0575 243  0.456  0.9778 
##  overweight - obese   0.000247 0.0575 243  0.004  0.9966 
## 
## area = Cerebellum:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.009050 0.0144 243  0.629  0.5299 
##  normal - obese       0.018253 0.0144 243  1.269  0.5299 
##  overweight - obese   0.009203 0.0144 243  0.640  0.5299 
## 
## area = CinguloOperc:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.008303 0.0470 243  0.177  0.9401 
##  normal - obese       0.004767 0.0470 243  0.101  0.9401 
##  overweight - obese  -0.003537 0.0470 243 -0.075  0.9401 
## 
## area = Default:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.022142 0.0308 243  0.720  0.5875 
##  normal - obese       0.038847 0.0308 243  1.263  0.5875 
##  overweight - obese   0.016706 0.0308 243  0.543  0.5875 
## 
## area = DorsalAttn:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.006997 0.0407 243  0.172  0.9966 
##  normal - obese       0.006825 0.0407 243  0.168  0.9966 
##  overweight - obese  -0.000172 0.0407 243 -0.004  0.9966 
## 
## area = FrontoParietal:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.002437 0.0332 243  0.073  0.9416 
##  normal - obese       0.030411 0.0332 243  0.915  0.6009 
##  overweight - obese   0.027974 0.0332 243  0.842  0.6009 
## 
## area = MedialParietal:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.012066 0.0814 243  0.148  0.8822 
##  normal - obese       0.039216 0.0814 243  0.482  0.8822 
##  overweight - obese   0.027149 0.0814 243  0.334  0.8822 
## 
## area = None:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight -0.004149 0.0814 243 -0.051  0.9594 
##  normal - obese      -0.026515 0.0814 243 -0.326  0.9594 
##  overweight - obese  -0.022367 0.0814 243 -0.275  0.9594 
## 
## area = Putamen:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.025571 0.0470 243  0.544  0.6680 
##  normal - obese       0.045746 0.0470 243  0.974  0.6680 
##  overweight - obese   0.020175 0.0470 243  0.429  0.6680 
## 
## area = Smhand:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.037030 0.0407 243  0.910  0.5455 
##  normal - obese       0.044031 0.0407 243  1.082  0.5455 
##  overweight - obese   0.007001 0.0407 243  0.172  0.8635 
## 
## area = Smmouth:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.027899 0.0575 243  0.485  0.7582 
##  normal - obese       0.045632 0.0575 243  0.793  0.7582 
##  overweight - obese   0.017733 0.0575 243  0.308  0.7582 
## 
## area = Thalamus:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.017180 0.0407 243  0.422  0.9464 
##  normal - obese       0.019916 0.0407 243  0.489  0.9464 
##  overweight - obese   0.002736 0.0407 243  0.067  0.9464 
## 
## area = VentralAttn:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight -0.007346 0.0575 243 -0.128  0.9523 
##  normal - obese      -0.010793 0.0575 243 -0.188  0.9523 
##  overweight - obese  -0.003447 0.0575 243 -0.060  0.9523 
## 
## area = VentralDiencephalon:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.023993 0.0575 243  0.417  0.9804 
##  normal - obese       0.022575 0.0575 243  0.392  0.9804 
##  overweight - obese  -0.001417 0.0575 243 -0.025  0.9804 
## 
## area = Visual:
##  contrast             estimate     SE  df t.ratio p.value
##  normal - overweight  0.012462 0.0210 243  0.593  0.8674 
##  normal - obese       0.000767 0.0210 243  0.037  0.9709 
##  overweight - obese  -0.011695 0.0210 243 -0.557  0.8674 
## 
## Degrees-of-freedom method: df.error 
## P value adjustment: fdr method for 3 tests
cent.gls <- gls(centrality ~ mods, data = warp)
cent.emm <- emmeans(cent.gls, "mods")
cent.fac <- update(cent.emm, levels = list(
                group = c("normal", "overweight", "obese"), 
                area = c("","Amygdala", "Auditory", "BrainStem", "Caudate", "Cerebellum", "CinguloOperc",
                         "Default", "DorsalAttn", "FrontoParietal", "MedialParietal", "None", "Putamen",
                         "Smhand", "Smmouth", "Thalamus", "VentralAttn", "VentralDiencephalon", "Visual")))

contrast(cent.fac, "pairwise", by = "area", adjust = "fdr")
## area = :
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -2.21e-04 0.002159 243 -0.103  0.9184 
##  normal - obese       3.09e-04 0.002159 243  0.143  0.9184 
##  overweight - obese   5.30e-04 0.002159 243  0.246  0.9184 
## 
## area = Amygdala:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight  1.06e-03 0.003053 243  0.348  0.8321 
##  normal - obese       1.71e-03 0.003053 243  0.561  0.8321 
##  overweight - obese   6.48e-04 0.003053 243  0.212  0.8321 
## 
## area = Auditory:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -5.20e-04 0.003053 243 -0.170  0.9926 
##  normal - obese      -5.48e-04 0.003053 243 -0.180  0.9926 
##  overweight - obese  -2.84e-05 0.003053 243 -0.009  0.9926 
## 
## area = BrainStem:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight  2.75e-05 0.001080 243  0.025  0.9921 
##  normal - obese       1.07e-05 0.001080 243  0.010  0.9921 
##  overweight - obese  -1.67e-05 0.001080 243 -0.016  0.9921 
## 
## area = Caudate:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -3.74e-04 0.002159 243 -0.173  0.9541 
##  normal - obese       1.24e-04 0.002159 243  0.058  0.9541 
##  overweight - obese   4.99e-04 0.002159 243  0.231  0.9541 
## 
## area = Cerebellum:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -1.45e-05 0.000540 243 -0.027  0.9786 
##  normal - obese      -1.40e-04 0.000540 243 -0.259  0.9786 
##  overweight - obese  -1.26e-04 0.000540 243 -0.233  0.9786 
## 
## area = CinguloOperc:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -5.63e-05 0.001763 243 -0.032  0.9746 
##  normal - obese       3.26e-04 0.001763 243  0.185  0.9746 
##  overweight - obese   3.82e-04 0.001763 243  0.217  0.9746 
## 
## area = Default:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -4.72e-04 0.001154 243 -0.409  0.8564 
##  normal - obese      -6.81e-04 0.001154 243 -0.590  0.8564 
##  overweight - obese  -2.09e-04 0.001154 243 -0.181  0.8564 
## 
## area = DorsalAttn:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -6.86e-05 0.001527 243 -0.045  0.9642 
##  normal - obese       1.06e-04 0.001527 243  0.069  0.9642 
##  overweight - obese   1.74e-04 0.001527 243  0.114  0.9642 
## 
## area = FrontoParietal:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -3.91e-05 0.001247 243 -0.031  0.9750 
##  normal - obese      -5.01e-04 0.001247 243 -0.402  0.9750 
##  overweight - obese  -4.62e-04 0.001247 243 -0.371  0.9750 
## 
## area = MedialParietal:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight  1.18e-04 0.003053 243  0.039  0.9692 
##  normal - obese      -6.35e-04 0.003053 243 -0.208  0.9692 
##  overweight - obese  -7.53e-04 0.003053 243 -0.247  0.9692 
## 
## area = None:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -5.95e-04 0.003053 243 -0.195  0.8456 
##  normal - obese       9.49e-04 0.003053 243  0.311  0.8456 
##  overweight - obese   1.54e-03 0.003053 243  0.506  0.8456 
## 
## area = Putamen:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight  2.65e-04 0.001763 243  0.150  0.9026 
##  normal - obese      -2.16e-04 0.001763 243 -0.123  0.9026 
##  overweight - obese  -4.81e-04 0.001763 243 -0.273  0.9026 
## 
## area = Smhand:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -5.56e-04 0.001527 243 -0.364  0.8934 
##  normal - obese      -7.60e-04 0.001527 243 -0.498  0.8934 
##  overweight - obese  -2.05e-04 0.001527 243 -0.134  0.8934 
## 
## area = Smmouth:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -2.54e-04 0.002159 243 -0.117  0.9066 
##  normal - obese      -1.00e-03 0.002159 243 -0.463  0.9066 
##  overweight - obese  -7.46e-04 0.002159 243 -0.346  0.9066 
## 
## area = Thalamus:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight  2.36e-04 0.001527 243  0.155  0.9909 
##  normal - obese       2.53e-04 0.001527 243  0.166  0.9909 
##  overweight - obese   1.74e-05 0.001527 243  0.011  0.9909 
## 
## area = VentralAttn:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight  1.11e-03 0.002159 243  0.515  0.8596 
##  normal - obese      -3.82e-04 0.002159 243 -0.177  0.8596 
##  overweight - obese  -1.49e-03 0.002159 243 -0.692  0.8596 
## 
## area = VentralDiencephalon:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight -4.67e-04 0.002159 243 -0.216  0.9709 
##  normal - obese      -7.89e-05 0.002159 243 -0.037  0.9709 
##  overweight - obese   3.88e-04 0.002159 243  0.180  0.9709 
## 
## area = Visual:
##  contrast             estimate       SE  df t.ratio p.value
##  normal - overweight  2.07e-05 0.000788 243  0.026  0.9791 
##  normal - obese       3.85e-04 0.000788 243  0.489  0.9660 
##  overweight - obese   3.65e-04 0.000788 243  0.463  0.9660 
## 
## Degrees-of-freedom method: df.error 
## P value adjustment: fdr method for 3 tests

Nice function to make a pretty table

library(mosaic)
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 3.6.2
## Loading required package: ggformula
## Loading required package: ggstance
## Warning: package 'ggstance' was built under R version 3.6.2
## 
## Attaching package: 'ggstance'
## The following objects are masked from 'package:ggplot2':
## 
##     geom_errorbarh, GeomErrorbarh
## 
## New to ggformula?  Try the tutorials: 
##  learnr::run_tutorial("introduction", package = "ggformula")
##  learnr::run_tutorial("refining", package = "ggformula")
## Loading required package: mosaicData
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Note: If you use the Matrix package, be sure to load it BEFORE loading mosaic.
## 
## Have you tried the ggformula package for your plots?
## 
## Attaching package: 'mosaic'
## The following object is masked from 'package:Matrix':
## 
##     mean
## The following objects are masked from 'package:psych':
## 
##     logit, rescale
## The following object is masked from 'package:plyr':
## 
##     count
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## The following object is masked from 'package:purrr':
## 
##     cross
## The following object is masked from 'package:ggplot2':
## 
##     stat
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum
tablr<-function(Y,x,D){
  Q<-favstats(Y ~ x, data = D)
  Q.stat <- Q[, c("x", "n", "mean", "sd")]
  colnames(Q.stat)<-c("test","n", "mean", "sd")
  a<-match.call()[2]
  return(Q.stat)
}

quickr<-function(X, Y, Z, X2, Y2, Z2, X3){
  a<-merge(X, Y, by="test")
  b<-merge(a, Z, by="test")

  c<-merge(X2, Y2,by="test")
  d<-merge(c, Z2, by="test")
  # e<-merge(d, X3, by="test")
  d$test<-revalue(d$test, c( "no" = "Total",
                             "ov" = "Total",
                             "ob" = "Total"))


  e<-rbind(b,d)
  f<-merge(e,X3, by="test")
  drops <- c("Y")
  f<-f[ , !(names(f) %in% drops)]
  # f<-rbind(e,X3)
  
  library("plyr")
  f$test<-revalue(f$test, c("Am. Indian/Alaskan Nat."=  "Native American",            
                                  "Asian/Nat. Hawaiian/Othr Pacific Is."=  "Asian, Native Hawaiian, or Pacific Islander",
                                  "Black or African Am."= "Black or African American",
                                  "More than one"= "More than one race",                      
                                  "Unknown or Not Reported"= "Unknown or chose not to report",
                                  "White"= "White",
                                  "Total"="Total"))
    f$test <- factor(f$test, levels = c( "Native American",            
                                  "Asian, Native Hawaiian, or Pacific Islander",
                                  "Black or African American",
                                  "More than one race",                      
                                  "Unknown or chose not to report",
                                  "White",
                                  "Total"))

  return(f)
}

sexr<-function(X,Y){
  a <- table(X,Y)
  b <- prop.table(a,margin=2)
  c<-round(b*100, 1)
  d<-data.frame(c)
  e<-subset(d, d$Y == "F")
  e<-rename(e, c("X"="test"))
  A<-table(Y)
  B <- prop.table(A)
  C<-round(B*100, 1)
  D<-data.frame(C)
  E<-subset(D, D$Y == "F")
  test <- rep("Total",length(1))
  G <- cbind(test, E)
  f<-rbind(e,G)
  # return(rename(f, c("X"="test")))
  return(f)

}

“Native American”, “Asian, Native Hawaiian, or Pacific Islander”,“Black or African American”,“More than one race”,“Unknown or chose not to report”,“White”,“Total” X Y Freq

demo_vars<-c("Race","BMI","HbA1C","Age_in_Yrs","group", "Gender")
dems<-dat[demo_vars]

groupSex<-table(dems$group, dems$Gender)
groupSex <- prop.table(groupSex,margin=2)
round(groupSex*100, 1)
##     
##         F    M
##   no 52.6 38.4
##   ov 24.6 39.0
##   ob 22.8 22.7
#dems<-na.omit(dems)
dem_no<-subset(dems, dems$group == "no")
dem_no$group<-factor(dem_no$group)

dem_ov<-subset(dems, dems$group == "ov")
dem_ov$group<-factor(dem_ov$group)


dem_ob<-subset(dems, dems$group == "ob")
dem_ob$group<-factor(dem_ob$group)

Making a nice table

## The following `from` values were not present in `x`: ov, ob
## The following `from` values were not present in `x`: no, ob
## The following `from` values were not present in `x`: no, ov
colz<-c("Race/Ethncity",
        "n","mean","SD",
        "n","mean","SD",
        "n","mean","SD",
        "%")
digz<-c(0, 2, 2, 0, 2, 2,0, 2, 2)

Make the nice table

## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
Table 1: Descriptive statistics by reported race and BMI group
BMI
HbA1C
Age at scan
Female
Race/Ethncity n mean SD n mean SD n mean SD %
Average BMI
Native American 0 NaN NA 0 NaN NA 0 NaN NA 0.0
Asian, Native Hawaiian, or Pacific Islander 19 21.48 2 12 5.16 1 19 26.53 4 12.2
Black or African American 23 22.73 1 12 5.40 0 23 28.22 4 16.7
More than one race 6 22.63 2 6 5.30 0 6 25.67 3 5.6
Unknown or chose not to report 0 NaN NA 0 NaN NA 0 NaN NA 0.0
White 108 22.18 2 73 5.18 0 108 28.81 4 65.6
Total 156 22.20 2 103 5.21 0 156 28.33 4 57.7
High BMI
Native American 1 29.23 NA 1 5.90 NA 1 35.00 NA 2.4
Asian, Native Hawaiian, or Pacific Islander 6 25.86 1 2 5.00 0 6 25.33 5 7.1
Black or African American 17 27.93 2 12 5.41 0 17 28.53 3 21.4
More than one race 4 26.31 1 2 5.20 0 4 22.75 0 4.8
Unknown or chose not to report 2 27.69 2 2 5.25 0 2 29.00 4 0.0
White 79 27.16 1 57 5.19 0 79 28.51 3 64.3
Total 109 27.20 2 76 5.23 0 109 28.19 4 38.5
Very high BMI
Native American 0 NaN NA 0 NaN NA 0 NaN NA 0.0
Asian, Native Hawaiian, or Pacific Islander 0 NaN NA 0 NaN NA 0 NaN NA 0.0
Black or African American 14 34.73 5 7 6.19 2 14 30.57 4 25.6
More than one race 2 32.05 2 2 5.70 0 2 30.00 0 5.1
Unknown or chose not to report 4 32.53 3 4 5.10 0 4 26.00 2 5.1
White 58 34.30 3 35 5.25 0 58 29.26 4 64.1
Total 78 34.23 4 48 5.39 1 78 29.35 4 50.0